Lanthanide-doped MoS2 with enhanced oxygen reduction activity and biperiodic chemical trends

Molybdenum disulfide has broad applications in catalysis, optoelectronics, and solid lubrication, where lanthanide (Ln) doping can be used to tune its physicochemical properties. The reduction of oxygen is an electrochemical process important in determining fuel cell efficiency, or a possible environmental-degradation mechanism for nanodevices and coatings consisting of Ln-doped MoS2. Here, by combining density-functional theory calculations and current-potential polarization curve simulations, we show that the dopant-induced high oxygen reduction activity at Ln-MoS2/water interfaces scales as a biperiodic function of Ln type. A defect-state pairing mechanism, which selectively stabilizes the hydroxyl and hydroperoxyl adsorbates on Ln-MoS2, is proposed for the activity enhancement, and the biperiodic chemical trend in activity is found originating from the similar trends in intraatomic 4f–5d6s orbital hybridization and interatomic Ln–S bonding. A generic orbital-chemistry mechanism is described for explaining the simultaneous biperiodic trends observed in many electronic, thermodynamic, and kinetic properties.


(A) Supplementary Structures and Stabilities of Ln-MoS2
The possible doping sites of Ln in MoS2 include Mo, S, and interstitial sites. After systematic structural calculations, we find that the Ln substituent at a Mo site is the most stable, and the derived Ln-MoS2 structure is stable against the structural relaxation. Apart from the stable Mo-site substitution, the other two dopant configurations are both unstable. The S-substituted Ln-MoS2 cannot resist against any surface adsorption during the structural relaxation, where the Ln dopant can be automatically extracted out of the lattice by a surface adsorbate ( Figure S1). During the structural relaxation of the interstitial-doped Ln-MoS2, the interstitial Ln atom even automatically destroy the local crystal structure.

(B) Supplementary Raman Spectra of Ln-MoS2
The Raman spectra of pristine MoS2 and Ln-MoS2 (e.g., Er-MoS2) have been measured in experiment [1]. It is necessary to study the effect of Ln dopant on the Raman spectra using accurate DFT calculation methods [2,3], which can not only further prove the rationality of the used structural model here, but also guide more detailed experimental characterizations in the future. In the calculated Raman spectrum for the pristine MoS2 (Figure S3), the positions of the two main Raman-active peaks, i.e., the in-plane E 1 2g mode at 384.5 cm -1 and the out-of-plane A1g mode at 406.7 cm -1 , are consistent with the experimental peaks at 384.2 ~ 385 and 403.8 ~ 404.9 cm -1 [4,5], respectively. After Er doping, the two main Raman-active peaks in the calculated spectrum are both red shifted down to 374.8 and 390.0 cm -1 (Figure S3), which is consistent with the experimental reports [1]. To explore the reason for such Raman-peak redshift caused by Er doping, we further calculate the Raman spectrum of the pristine MoS2 with the lattice constant expanded to be the same as that of Er-MoS2. The results show that the lattice expansion caused by Er dopant only has a partial contribution to the redshift, and the chemical effect of Er (e.g., changing the bond strength) should have an additional contribution. Furthermore, there are some vibraional hybridizations [6] occuring in the Raman-active modes, where each mode contains more or less both of the in-plane and outof-plane vibrations. It is such vibrational hybridization that has made the higher Raman peak in pristine MoS2 separated into many much lower peaks in both Er-MoS2 and expanded MoS2. The La and Lu dopants bring the largest and smallest lattice expansions to the MoS2 lattice, respectively, and the Raman spectra of pristine MoS2, Er-MoS2, La-MoS2, and Lu-MoS2 are compared in Figure  S4, where the qualitative correspondence between redshift of the lower peak and the excess area can be clearly seen.

(C) Supplementary Electronic Properties of Ln Elements in Various States
When in the state of atomic gas, the ground state valence electron configurations of La, Ce, Gd, and Lu atoms are 4f n-3 5d 1 6s 2 (n is the total number of 4f, 5d, and 6s electrons), and those of other Ln atoms are 4f n-2 6s 2 . The sublimation heats of elemental Ln metals [7] exhibit a biperiodic chemical trend with respect to Ln type (Figure S5a), which indicates the same trend for the interatomic bond strength in Ln metals. Similar biperiodic trend is also observed in the third and fourth ionization potentials of Ln atoms (IP, Figure S5b) [7], both of which can reflect the ionization capacity of 4f electrons. We make a weighted linear summation of IP3 + 0.2×IP4, where the coefficients are first obtained by roughly fitting to the formation energies of Ln dopants, and then optimized to be 0.2 for a good trend presentation. The overall biperiodic trend of IP3 + 0.2×IP4 is compared with the variation trend in the formation energy of Ln-MoS2 (see main text, Figure   1c). Furthermore, the biperiodic trend also appears in the homolytic LnF bond energies of LnF3 molecules (1/3LnF3 → 1/3Ln + F, Figure S6) [8]. It is the common important role of 4f orbitals that have resulted in such coincident biperiodic trends in various physical properties of Ln in different states (e.g., atom, metal, and compound molecule).  [7]. We also make a weighted linear summation of IP (IP3 + 0.2×IP4) to compare with the variation trend in the formation energy of Ln-MoS2 (see main text, Figure 1c). Source data are provided as a Source Data file. Figure S6. The homolytic LnF bond energies of LnF3 molecules. The reaction equation is 1/3LnF3 → 1/3Ln + F [8]. The data for La, Ce, Pr, Nd, Eu, Gd, Dy, Ho, Er, and Tm are deduced from experimental measurements, and those for Pm, Sm, Tb, Yb, and Lu are obtained from quantum-chemistry calculations. Source data are provided as a Source Data file.
To clearly understand the intrinsic characters of atomic orbitals here, we perform the fullpotential all-electron calculations for Ln atoms (residing at the center of a sphere with radius of 100 Å) using the pseudopotential generator (i.e., the LD1 module) of the Quantum ESPRESSO code package [9]. The obtained energy levels and radial wavefunctions (rϕ(r), r is the radial coordinate) of the 4f, 5d, and 6s orbitals are shown in Figure S7a and S7b. It can be seen that the 4f orbital of a Ln atom is clearly lower in energy and much localized in spatial distribution than those of the 5d and 6s orbitals. Thus, in a solid, it is the 5d and 6s orbitals on a Ln atom that can sufficiently bond with the orbitals of surrounding atoms. However, the highly-localized 4f orbitals have relatively weak bonding capability and mainly contribute to the magnetic moments of Ln-MoS2 ( Figure S7c), which also present a clear Hund-rule trend. The energy cost for the transition of 4f electrons up to the 5d and 6s orbitals [10] exhibit a biperiodic chemical trends with respect to the Ln type ( Figure S7d). This can be well explained by an atomic-orbital mechanism: (1) From La to Gd, the occupation of 4f orbitals approaches the half-filling state, the strong exchange attraction between electrons with the same spin orientation (i.e., Hund rule) and the increased ion-electron attraction both contribute to the increased electronic-transition energy; (2) When crossing the half-filling boundary from Eu to Gd, the electronic-transition energy decreased abruptly due to the additional 4f electron has no exchange attraction with other seven 4f electrons with different spin orientations; (3) From Gd to Yb, the Hund rule is obeyed again by the added 4f electrons, and the ground-state is more and more stabilized by the increased electronic exchange attraction, as well as the increased ion-electron attraction.

S-8
The biperiodic chemical trend in the energy cost for the transition of 4f electrons up to 5d and 6s indicates that there should be a similar trend in the hybridization capability of 4f with 5d and 6s, which makes part of 4f electrons excited on to the later two orbitals. The exothermic interatomic bonding tend to increase with the degree of such endothermic intraatomic orbital hybridization [11]. It is such intrinsic orbital mechanism that has led to the general biperiodic trends ubiquitously appearing in the sublimation heats of Ln metals (Figure S5a), homolytic bond energies of LnF3 molecules (Figure S6), and formation energies of Ln dopants in MoS2 (see main text, Figure 1c).

(D) Supplementary Details of ORR Steps
As the initial step of an ORR process, O2 molecule first diffuses from the bulk electrolyte to the material-electrolyte interface, and then be adsorbed on an active S atom close to the Ln dopant.
This process can be expressed as where O2(aq), O2(dl), and O2 * represent the O2 molecule in aqueous solution, electrical double layer near the surface, and adsorbed state, respectively. ki and k-i are the rate constants for the forward and reverse reaction steps.
After being adsorbed on the Ln-MoS2 surface, O2 can be reduced into H2O by electrochemical protonation steps, for which there are two possible mechanisms: where *A and *B represent two adjacent surface S atoms closest to the Ln dopant. For OOH * , in addition to being electrochemically reduced to H2O in Eq. (4) above, it can also be chemically dissociated along the path of (10) OOH * A + *B S-12

(E) Supplementary Gibbs Free Energy of Reaction
The Gibbs free energy of reaction for each elemental ORR step (∆G) is calculated using the computational hydrogen electrode (CHE) model [12]: where ∆ε0 is the change in electronic energy; ∆εzpe is the change in vibrational zero-point energy; T is temperature (= 298.15 K here); and ∆s is the entropy change. The last term is the contribution of electrode potential (U, referring to the reversible hydrogen electrode) and only needs to consider in the steps involving electron transfer. The entropies of clean/adsorbed surfaces are derived from the calculated vibrational frequencies, and those of gaseous molecules are collected from experimental data [13]. Liquid H2O and gaseous H2 are used as reference states, and the free energy of liquid H2O is derived from that of the gaseous H2O at 0.032 bar (i.e., the saturated-vapor-pressure of water at 298.15 K). The free energy of O2 is obtained from the total reaction O2(g) + 2H2(g) → 2H2O(l) with the reaction free energy of -4.92 eV at 1 bar and 298.15 K.

(F) Supplementary Definition of Adsorption Free Energies
The adsorption free energies (

(H) Supplementary Energies and Electronic Structures for Adsorbed Ln-MoS2
The ∆Gadss for different adsorbates on Ln-MoS2 surfaces, except for ∆Gads(O2), clearly exhibit the biperiodic trends ( Figure S18).   The adsorption structures and the differential electronic densities induced by the adsorption of ORR intermediates on the pristine MoS2 and Ln-MoS2 surfaces are shown in Figure S20. It can be seen that the length of SO bond in the OH * system is considerably shortened from 1.81 down to 1.65 Å after Ln doping, while that in the O * system has been barely changed. In addition, the adsorption state of OOH has a special change from physical adsorption to chemical adsorption after Ln doping. These changes in bond lengths are in accordance with the selectively enhanced adsorption of OH and OOH by Ln doping, as discussed in the main text. From the differential electronic densities on Ln-MoS2, we can clearly observe the highly hybridized atomic orbitals around the SO bond in O * , OH * , and OOH * on Ln-MoS2, with a higher hybridization degree than that on pristine MoS2, indicating the stronger covalent SO bonding on Ln-MoS2.  The ΔGads(OH) of Pt (111) is 0.80 eV [12]. Source data are provided as a Source Data file.   Taking Sm-MoS2 as a representative, the activation energy barriers (Ea,i) of all ORR steps are calculated with different water structures using the CI-NEB method [14], for which the results are shown in the main text ( Figure 4a) and the atomic structures for different reaction paths are shown in Figure S30. It can be seen that the dissociation energy barrier of O2 * is as high as 1.26 ~ 1.50 eV and then is difficult to overcome at room temperature, while the transition from O2 * to OOH * only requires a small energy barrier of 0.05~0.36 eV. Similarly, the Ea,10 of OOH * dissociation (0.42~0.68 eV) is also considerably larger than that (0.01~0.03 eV) of the associative transition from OOH * to O * . Therefore, the ORR process on Ln-MoS2 surface will take the OOH association pathway. Inspired by Hansen's work [15], the microkinetic model for the ORR process on Ln-MoS2 surface can be constructed. The rate equations with the involved species (Eq. (1) ~ (6) in Section D) are given below: x O 2 (aq) = 2.34 × 10 -5 , where θ is the coverage of ORR intermediates, t is time, and x O 2 (aq) and x O 2 (dl) are the mole fraction of O2 in aqueous solution and electric double layer, respectively. These rate equations are solved at steady state, and then the turn over frequency and current density are finally calculated.
The rate constants k1 and k-1 of the diffusion of O2 molecule through the Nernstian diffusion layer on the rotating disk electrode (RDE) are calculated using [15] k 1 = k -1 = According to the transition-state theory, k2 for the adsorption process of O2 is calculated by where k B is Boltzmann constant, T is Temperature, h is Planck constant, and G a,2 is the activation free energy. Due to the weak adsorption of O2 on Ln-MoS2, G a,2 is manly contributed by the reorientation effect of solvent molecules, which is calculated to be 0.07 ~ 0.17 eV (see Figure S16 and Figure 2c in the main text) here. It has a higher magnitude (0.22 ~ 0.28 eV) on Pt [15] due to the stronger water adhesion. In our testing calculations using G a,2 from 0.07 to 0.28 eV, Uhalf has negligible variation (<0.002 V, Figure S31), because O2(dl) → O2 * is not the ORR-limiting step (see main text, Figure 4e). The ki of the i'th electrochemical step is calculated by where G a,i 0 is the activation free energy at the equilibrium potential (U i 0 ) of the i'th step, α i is the transfer coefficient (set to be 0.5), and U is the electrode potential. U i 0 is calculated by where ΔG i 0 is the free-energy change of the i'th reaction step at 0 V. The reverse-reaction rate constant k -i is calculated by exploiting the equilibrium constant (K i ): The current density (j) is calculated by where 4 is the number of transferred electrons, ρ is the surface density of active sites, and TOF O 2 is the turnover frequency of O2. The kinetic parameters with constant values are summarized in Table S1.   [15] for the O2 *  OOH * step is shown in the main text (Figure 4b), which can be used to estimate the Ea,3's of other Ln-MoS2 systems when having calculated the corresponding reaction energies (ΔE3). Due to the important role of this step in determine the polarization curves (see Figure 4e in the main text), we further investigate the structural mechanism underlying the Ea,3 with the example of Sm-MoS2. In Figure S24 and S32a, we can see that different water configurations will result in different distances between O2 * and its nearest H atom (dO-H), which Ea,3 closely correlates with, and the average Ea,3 is 0.19 eV. We further perform an AIMD simulation (300 K, 40 ps, 0.5 fs/step) on the structure of O2 * with a H atom nearby, the variation of dO-H ( Figure S32b) clearly shows that apart from freely moving within the water layer in most of the time, the H neighbor frequently gets close to the O2 molecule (1.8~2.0 Å), which can facilitate the following protonation reaction of O2 * . Thus, the Sm-MoS2/water interface with a dO-H around 2 Å (e.g., 1.86 Å) should be chosen to derive an appropriate Ea,3, and a too short (e.g., 1.6 Å) or too long (e.g., 2.7 Å) will result in an underestimated or overestimated reaction barrier ( Figure S32a). Furthermore, although the adsorption of O2 * on Sm-MoS2 surface is relatively weak, it does not escape from the Sm-MoS2/water interface during the AIMD simulation, indicating its preferred stability to act as an initial surface reactant for the whole highly-active ORR process. The variation of dOH in the AIMD simulation at 300 K for Sm-MoS2/water interface with an O2 * on surface and an additional H atom in water. Source data are provided as a Source Data file.

S-38
For the OOH *  O * step, the Ea,4 of Sm-MoS2 is lower than 0.03 eV as shown in the main text ( Figure 4a) and the U 4 0 is around 2.72 V, indicating that this step will proceed spontaneously when the electrode potential is lower than 1.23 V. Thus, we can safely assume that the Ea,4's for other Ln-MoS2 surfaces all equal 0 eV. For the O *  OH * step, the calculated Ea,5's of La-, Ce-, Pr-, Sm-, Dy-, and Yb-MoS2 are very close to each other around 0.12 eV (Figure S33a), except for Lu-MoS2 (~ 0.09 eV), thus 0.12 eV is adopted for other Ln-MoS2 surfaces. The Brønsted-Evans-Polanyi relationship [15] for the Ea,6 (OH *  H2O) values is also shown in Figure S33b. Although there is a certain deviation of the as-calculated data from this linear relationship, the total ORR reaction rate has a very small dependence on Ea,6 (i.e., having a XRC,6 close to 0, see Figure 4e in the main text). Therefore, a variation of 0.1 ~ 0.4 eV in Ea, 6 will cause negligible change in the reaction rate. The simulated current-potential polarization curves for all the Ln-MoS2 systems are shown in Figure S34.  The effects of RDE rotation speed (within a typical range from 900 to 3200 rpm) and temperature (within a typical range from 25 to 60 ℃) on the simulated polarization curves are evaluated here (Figure S37), where the representative Sm-MoS2 system is used. Both the rotation speed and temperature will affect the current density mainly through influencing the diffusion of O2 from the bulk electrolyte to the material-electrolyte interface. It can also be seen from Figure   S37 that although they can bring obvious change in the absolute value of current density at surface potential of ≤ 0.9 V, they only slightly shift the Uhalf by ~ 0.01 V. The above observations on the effects of RDE rotation speed and temperature actually are quite consistent in essence with the kinetic rate control analysis in the main text (Figure 4e). The effect of temperature can be further decomposed into reducing the initial concentration of O2 in water (x O 2 (aq), from 2.45×10 -5 at 25 ℃ to 1.58×10 -5 at 60 ℃ [18]) and increasing the reaction rate constant ( Figure S38). The temperature dependence of k1 and k-1 mainly originate from two factors (also see the formula for microkinetic model above): (1) the increase of O2 diffusion coefficient D O 2 from 2.4×10 -5 at 25 ℃ to 4.6×10 -5 cm 2 ·s -1 at 60 ℃, and (2) the decrease in water kinematic viscosity ν from 1.01×10 -6 at 25 ℃ to 0.47×10 -6 m 2 ·s -1 at 60 ℃ [18]. Apart from k1 and k-1, the reaction rates for the transitions between ORR intermediates at different temperatures are derived from the transition-state theory (see the formula for ki above). In Figure S38, the dependence of current density on the individual thermal changes of rate constant and O2 concentration is quantitatively disentangled. In addition, the effect of RDE rotation speed is attributed to the changed Nernstian diffusion layer and the associated changes in k1 (see its formula above). Figure S38. The thermal-effect mechanism disentangling for the polarization curve of the representative Sm-MoS2 system. The temperatures of 40 and 60 ℃ are considered, and the curves denoted with keys of "40 ℃" ("60 ℃") , "40-k" ("60-k"), and "40-x O 2 (aq) " ("60-x O 2 (aq) ") correspond to the simulated overall curves, the curves simulated only with the rate constants changed in accordance with temperature, and the curves simulated only with x O 2 (aq) changed, respectively. Source data are provided as a Source Data file.

(L) Supplementary Numerical Tests Validating the DFT Accuracy
When dealing with Ln elements that have highly localized 4f orbitals, there is the electronic self-interaction problem that we need to concern for the conventional PBE-GGA functional [30].
We use the efficient DFT plus Hubbard U method [31] to probe the influence of this problem on the accuracy of calculated surface reactivity. The tested adsorption free energies (ΔGads) of O and OH adsorbates on the active S sites of Ce-MoS2 and Gd-MoS2, which respectively have the minimum and maximum magnetic moments (see Figure S7c), are found nearly independent on the Ueff value (energy variation <0.03 eV, see Table S4). Furthermore, the effect of spin-orbit coupling (SOC) for the 4f orbitals is also tested, and its influence on the ΔGadss of adsorbates on active S site is also negligible (energy variation <0.01 eV, see Table S5). It is an indirect way for the Ln dopants to affect the ORR activity of MoS2 by changing the Ln-S bond character, and the 4f orbitals on Ln atom will hybridize with the higher-energy orbitals (e.g., 5d, 6s, and 6p), which contribute to the bond character and strength. The convergence testing calculations for ΔGadss using different cutoff energies and reciprocal-point mesh sizes are shown in Table S6 and S7, which clearly show that the settings used in this work (450 eV and 2×2×1) are stringent enough to yield accurate ΔGadss (energy convergence < 0.01 eV).  Regarding the structural model for Ln-MoS2, it is the monomer-Ln dopants observed in experimental samples (see Figure 1b in main text), where the inter-dopant distances are ~20 Å with a good dispersion. Therefore, such monomer dopant model is considered here for the study of ORR activity. To reveal the effect of Ln concentration on surface reactivity, the supercell-size effect on ΔGads is calculated (Table S8), which shows that an energy convergence of <0.02 eV can be achieved by the supercell of 4×2√3 used in this work. The weak supercell-size effect on Ln dopant can be also well reflected by the supercell-size dependence of Ef (Figure S39 and S40), and the supercell used in this work can guarantee an energy convergence of 0.1~0.2 eV (i.e., < 3%), and the biperiodic chemical trend in Ef has no obvious quantitative change when using a larger supercell (e.g., 6×2√3, Figure S40). Furthermore, the very weak magnetic coupling strength between two This situation is different from that of transition-metal doped nitrogenated graphene [35,36], where the active sites exactly reside on the transition-metal dopants and their magnetic coupling effect may play a role in the ORR activity. Therefore, according to the above comprehensive testing calculations for supercell-size effects, the structure model for Ln-MoS2 as used in this work (i.e., a Ln dopant in the 4×2√3 supercell) can well represent the realistic Ln-MoS2 samples with a high numerical accuracy.